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Abstract 

On the basis of statistical mechanics of the Q-Ising model, we formulate the Bayesian 
inference to the problem of inverse halftoning, which is the inverse process of repre- 
senting gray-scales in images by means of black and white dots. Using Monte Carlo 
simulations, we investigate statistical properties of the inverse process, especially, we 
reveal the condition of the Bayes-optimal solution for which the mean-square error 
takes its minimum. The numerical result is qualitatively confirmed by analysis of 
the infinite-range model. As demonstrations of our approach, we apply the method 
to retrieve a grayscale image, such as standard image Lenna, from the halftoned ver- 
sion. We find that the Bayes-optimal solution gives a fine restored grayscale image 
which is very close to the original. 
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1 Introduction 



In recent two or three decades, a considerable number of researchers have 
investigated various problems in information sciences, such as image restora- 
tion and error-correcting codes on the basis of the analogy between statistical 
mechanics and probabilistic information processing [1] . Especially, a lot of re- 
searchers have investigated various problems in image processing based on the 
Markov random fields [2,3,4,5]. In the field of the print technologies, many 
techniques of information processing have also developed. Particularly, the 
digital halftoning [7,8,9,10,11] is regarded as a key processing to convert a 
digital grayscale image to black and white dots which represents the original 
grayscale levels appropriately. On the other hand, the inverse process of the 
digital halftoning is referred to as inverse halftoning. The inverse halftoning 
is also important for us to make scanner machines to retrieve the original 
grayscale image by making use of much less informative materials, such as the 
halftoned binary dots. The inverse halftoning is 'ill-posed' in the sense that 
one lacks information to retrieve the original image because the material one 
can utilize is just only the halftoned black and white binary dots instead of 
the grayscale one. To overcome this difficulty, we usually introduce the 'regu- 
larization term' which compensates the lack of the information and regard the 
inverse problem as a combinatorial optimization [12,13]. Then, the optimiza- 
tion is achieved to find the lowest energy state via, for example, simulated 
annealing [14,15]. 

Besides the standard regularization theory, we can use the Bayesian approach. 
Under the direction of this approach, Stevenson [16] attempted to apply the 
maximum of a Posteriori (MAP for short) estimation to the problem of in- 
verse halftoning for a given halftone binary dots obtained by the threshold 
mask and the so-called error diffusion methods. However, there is few theoret- 
ical approach to deal with the inverse-halftoning from the view point of the 
Bayesian inference and statistical mechanics of information. 

In this study, on the basis of statistical mechanics of the Q-Ising model [17], we 
formulate the problem of inverse halftoning to estimate the original grayscale 
levels by using the information about both the halftoned binary dots and 
the threshold mask. We reconstruct the original grayscale revels from a given 
halftoned binary image and the threshold mask so as to maximize the posterior 
marginal probability. Using Monte Carlo simulations, we investigate statistical 
properties of the inverse process, especially, we reveal the condition of the 
Bayes-optimal solution for which the mean-square error takes its minimum. 
The result of the simulation is supported by the analysis of the infinite-range 
model. In order to investigate to what extent the Bayesian approach is effective 
for realistic images, we apply the method to retrieve the grayscale levels of 
the 256-levels standard image henna from the binary dots. We find that the 
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Bayes-optimal solution gives a fine restored grayscale image which is very close 
to the original one. 

The contents of this paper are organized as follows. In the next section, we 
formulate the problem of inverse halftoning. We mention the relationship be- 
tween statistical mechanics of the Q-Ising model and Bayesian inference of the 
inverse halftoning. In the following section, we investigate statistical proper- 
ties of the Bayesian inverse halftoning by Monte Carlo simulations. Analysis of 
the infinite-range model supports the result of the simulations. We also show 
that the Bayes-optimal inverse halftoning is useful even for realistic images, 
such as the 256-level standard image henna. Last section is summary. 



2 The model 



We first define the model system to investigate the statistical performance 
of the Bayesian inference for the problem of inverse halftoning. As original 
grayscale images, which are converted to the black and white binary dots, 
we consider snapshots from a Gibbs distribution of the ferromagnetic Q-Ising 
model having the spin variables {£} = {Cx,y — 0, • • • , Q — l\x, y — 0, • • • , L — 
1}. Then, each image {£} being specified by the Hamiltonian = 
^En.n.(Ci,t/ — £x',y') 2 f°ll° ws the Gibbs distribution 



Pr({£}) = — exp 





1 




= ^exp 



(1) 



at temperature T s , where Z s is the partition function of the system and the 
summation X) n .n.(' ' ') runs over the sets of the nearest neighboring pixels lo- 
cated on the square lattice in two dimension. The ratio of strength of spin-pair 
interaction J s and temperature T s , namely, J s /T s controls the smoothness of 
our original image {£}. In Fig. 1 (left), we show a typical example of the snap- 
shots from the distribution (1) for the case of Q = 4, J s = 1 and T s = 0.5. 
The right panel of the Fig. 1 shows the 256-levels grayscale standard image 
Lenna with 400 x 400 pixels. We shall use the standard image to check the 
efficiency of our approach in the last part of this paper. 

In order to convert original grayscale images to the the black and white binary 
dots, we make use of the threshold array {M}. Each component M^i of the 
array {M} takes a non-overlapping integer and these numbers are arranged 
on the h m x h m squares as shown in Fig. 2 for h m = 2 (left) and for h m = 4 
(right). For general case of L m , we define the array as 
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Fig. 1. An original image as a snapshot from the Gibbs distribution of (1) having 
100 x 100 pixels for the case of Q = 4 (left). We set T s = 0.5, J = 1. The right panel 
shows a 256-levels standard image henna with 400 x 400 pixels. 






2 


3 


1 






8 


2 


10 


12 


4 


14 


6 


3 


11 


1 


9 


15 


7 


13 


5 



Fig. 2. The Bayer-type threshold arrays for the dither method with 2x2 (left) and 
with 4x4 (right). 



{M} = {M k>l = 0, 



Q-l 2{Q-l) 



LI 



1' LI 



1 ' 



k, I — 0, 1, • • • , L r 



(2) 



We should keep in mind that the definition (2) is reduced to {M} = {Mk y i = 
0, 1, • • • , Q — l\k, I — 0, 1, • • • , \[Q — 1} and the domain of each component of 
the threshold array becomes the same as that of the original image {£} for 
L 2 m = Q- 

In order to achieve a pixel-to-pixel map between each element of the threshold 
array, M XjV and the corresponding original grayscale pixel £ XjJ/ , we spread a lots 
of threshold arrays over the original image so as not to overlap any threshold 
array with one another. Then, we transform each original pixel C, X)V into the 
binary dot T x>y by 



T x,y — \£x,y M Xj y) . 



(3) 



Here we denned M x ^ y as the threshold value corresponding to the (x, y)-th pixel 
and 9{- • •) denotes the unit-step function. Halftone images generated by the 
dither method via (3) are shown in Fig. 3. We find that the left panel obtained 
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by the uniform threshold mask M x>y = 2 (V X)J/ ) is hard to be recognized as a 
grayscale image, whereas, the center panel obtained by the 2x2 Bayer-type 
threshold array might be recognized as just like an original image through our 
human vision systems (due to a kind of optical illusion). 

Obviously, the inverse process of the above halftoning is regarded as an ill- 
posed problem. This is because from (3), one can not determine the origi- 
nal image £ X;J/ (V X;J/ ) completely from a given set of t XjV (V XjJ/ ) and M XjV (V XiJ/ ). 
Then, the standard regularization theory [12,13] provides us a realistic break- 
through. In the theory, we introduce the so-called 'regularization term' that 
compensates the lack of the information to retrieve the original image. Then, 
we construct the energy function to be minimized to find the original image 
as the lowest energy state. For instance, Some recent progress based on the 
standard regularization theory is found in our paper [18]. 

The standard regularization theory is itself a general and powerful approach, 
nevertheless, we here use an alternative, namely, the Bayesian approach to 
solve the inverse problem. 




Fig. 3. The left panel shows a halftone image converted by the dither method using 
the uniform threshold M = 2 from the snapshot from a Gibbs distribution of the 
Q = 4 Ising model shown in Fig. 1 (left). The center panel shows a halftone image 
obtained by the dither method using the 2x2 Bayer-type threshold array from 
the same snapshot. The right panel shows a halftone image converted by the dither 
method using the 4x4 Bayer-type threshold array from the 256-level standard 
image henna with 400 x 400 pixels shown in Fig. 1 (right). 



In the Bayesian inverse digital halftoning, we attempt to restore the origi- 
nal grayscale image from a given halftone image by means of the so-called 
maximizer of posterior marginal (MPM for short) estimate. Then, we define 
{z} = {z XtV — 0, • • • , Q — l\x, y — 0, • • • , L — 1} as an estimate of the original 
image {£ } arranged on the square lattice and reconstruct the grayscale image 
on the bases of maximizing the following posterior marginal probability: 
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z x , y = arg max V Pr({z}|{r}) = argmaxPr (^ i2/ |{r}) , (4) 

ZX ' V ! \-L ZX ' V 

\ z JT :z x,y 

where the summation J2 Zx y ^{z}(' ' *) runs over an pixels except for the (x, y)-th 
and the posterior probability -P({-z}|{t}) is given by the Bayes formula: 

pr (u\\w\) - YiMl^MM} (5) 

Pr(MIW) "E W Pr({4)Pr({r}|{,}) ^ 



In this study, following Stevenson [16], we assume that the likelihood might 
have the same form as the halftone process of the dither method, namely, 

P ({r}|{£}) = U {XtV) S (r x , y , 6 (z x , y - M x , y )) , (6) 



where 5 (a, b) denotes a Kronecker delta and we should notice that the infor- 
mation on the threshold array {M} is available in addition to the halftone 
image {r}. Then, we choose the model of the true prior as 



Pr({^}) = — — exp 



-m n.n. 



(7) 



where Z m is a normalization factor. J and T are the so-called hyper-parameters. 
It should be noted that one can construct the Bayes-optimal solution if we 
assume that the model prior has the same form as the true prior, namely, 
J = J s and T rn = T s (what we call, Nishimori line in the research field of spin 
glasses [1]). 

From the viewpoint of statistical mechanics, the posterior probability Pr({^}|{r}) 
generates the equilibrium states of the ferromagnetic Q-Ising model whose 
Hamiltonian is given by 

H ({z}) = Jj2( z x,y~ z x',y') 2 , (8) 



under the constraints 

^x,y Tx,y = & { z x,y ~ ^x,y) ■ (9) 



Obviously, the number of possible spin configurations that satisfy the above 
constraints (9) is evaluated as Y[(x, y ) \Q T x, y — M x>y \ and this quantity is expo- 
nential order such as ~ a L2 (a: a positive constant). Therefore, the solution 
{z} to satisfy the constraints (9) is not unique and this fact makes the problem 
very hard. To reduce the difficulties, we consider the equilibrium state gen- 
erated by a Gibbs distribution of the ferromagnetic Q-Ising model with the 
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constraints (9) and increase the parameter J gradually from J = 0. Then, we 
naturally expect that the system stabilizes the ferromagnetic Q-Ising config- 
urations due to a kind of the regularization term (8). Thus, we might choose 
the best possible solution among a lot of candidates satisfying (9). 

/From the view point of statistical mechanics, the MPM estimate is rewritten 
by 

Zx,v = 9q{{Zx,v)), {Zx,y)=^2z Xiy PT({z}\{T}) (10) 

2 

where 0q(- • •) is the Q-generalized step function defined by 

e Q (x) = £ k{o (x - (k - £)) - e [x - (k + £)) }. (ii) 

Obviously, (z x y ) is a local magnetization of the system described by (8) under 
(9). 

2.1 Average case performance measure 

To investigate the performance of the inverse halftoning, we evaluate the mean 
square error which represents the pixel-wise similarity between the original and 
restored images. Especially, we evaluate the average case performance of the 
inverse halftoning through the following averaged mean square error 

^ = 7^E Pr ({0)E(^-^) 2 . (12) 

We should keep in mind that the a gives zero if all restored images are exactly 
the same as the corresponding original images. 



3 Results 

In this section, we first investigate the statistical properties of our approach 
to the inverse halftoning for a set of snapshots from a Gibbs distribution of 
the ferromagnetic Q-Ising model via computer simulations. We next analyti- 
cally evaluate the performance for the infinite-range model. Finally, we check 
the usefulness of our approach for the realistic images, namely, the 256-levels 
standard image henna. 
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3.1 Monte Carlo simulation 



We first carry out Monte Carlo simulations for a set of halftone images, which 
are obtained from the snapshots from a Gibbs distribution of the ferromagnetic 
Q = 4 Ising model with 100 x 100 pixels by the uniform threshold M x>y = 
2 (V XjJ/ ) and the 2x2 Bayer-type threshold arrays as shown in Fig. 2. In order 
to clarify the statistical performance of our method, we reveal the hyper- 
parameters J and T m dependence of the averaged mean square error a. We 
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Fig. 4. The mean square error as a function of T m . The original image is a snapshot 
from a Gibbs distribution of the Q = 4 ferromagnetic Ising model with 100 x 100 
pixels and T s = 1.0, J s = 1 and J = 1. The halftone images are obtained by the 
uniform and 2x2 Bayer- type arrays. 

plot the results in Fig. 4. These figures show that the present method achieves 
the best possible performance under the Bayes-optimal condition, that is, 
J = J s and T m = T s . We also find from Fig. 4 (the lower panel) that the limit 
T m — > oo leading up to the MAP estimate gives almost the same performance 
as the Bayes-optimal MPM estimate. 

This fact means that it is not necessary for us to take the T m — > limit 
when we carry out the inverse halftoning via simulated annealing. From the 
restored image in Fig. 5 (center), it is actually confirmed that the present 
method effectively works for the snapshot of the ferromagnetic Q-Ising model. 

It should be noted that the mean square error evaluated for the 2x2 Bayer- 
type array is larger than that for the M = 2 uniform threshold. This result 
seems to be somewhat counter-intuitive because the halftone image shown 
in the center panel of Fig. 3 seems to be much closer to the original image, 
in other words, is much informative to retrieve the original image than the 
halftone image shown in the left panel of the same figure. However, it could 
be understood as follows. The shape of each 'cluster' appearing in the original 
image (see the left panel of Fig. 1) remains in the halftone version (the left 
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Fig. 5. The left panel shows a Q = 4 grayscale image restored by the MPM estimate 
from the halftone image shown in Fig. 3 (left). The center panel shows a Q = 4 
grayscale image restored by the MPM estimate from the halftone image shown in 
Fig. 3 (center). The right panel shows a Q = 256 grayscale image restored by the 
MPM estimate from the halftone image shown in Fig. 3 (right). 

panel of Fig. 3), whereas, in the halftone image (the center panel of Fig. 3), 
such structure is destroyed by the halftoning process via the 2x2 Bayer- 
type array. As we found, in a snapshot of the ferromagnetic Q-Ising model 
at the inverse temperature J s /T s = 1, the large size clusters are much more 
dominant components than the small isolated pixels. Therefore, the averaged 
mean square error is sensitive to the change of the cluster size or the shape, and 
if we use the constant threshold mask to create the halftone image, the shape 
of the cluster does not change, whereas the high-frequency components vanish. 
These properties are desirable for us to suppress the increase of the averaged 
mean square error. This fact implies us that the averaged mean square error 
for the 2x2 Bayer-type is larger than that for the constant mask array and 
the performance is much worse than expected. 



3.2 Analysis of the infinite- range model 

In this subsection, we check the validity of our Monte Carlo simulations, 
namely, we analytically evaluate the statistical performance of the present 
method for a given set of the snapshots from a Gibbs distribution of the ferro- 
magnetic Q-Ising model in which each spin variable is located on the vertices 
of the complete graph. For simplicity, we first transform the index from (x, y) 
to i so as to satisfy i — x + Ly + 1. Then, the new index i runs from i = 1 
to L 2 - 1 = N. For this new index of each spin variable, we consider the 
infinite-range version of true prior and the model as 



Pr({£}) 



g 2N ^2i<j(& €j) 2 



z s 



Pr({z}) 




(13) 
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where the scaling factors 1/N appearing in front of the sums Ei<j(' ' ') are 
needed to take a proper thermodynamic limit. We also set /3 S = J s /T s and 
f3 m = J/T m for simplicity. Obviously, the thermodynamics of the system {£} 
is determined by the following magnetization: 



1 ^ Egroeexp[2^m e-^ 2 ] 
m ° " ^ S exp[2/3 s m £ - AC 2 ] ' 

On the other hand, the magnetization for the system {z} having disorders {t} 
and {£} is given explicitly as 



1 



171 : 



T Q-1 ( Eta z ^ mmZ -" mz2s ( e (S-M),e(z-M)) \ ofrmot-fl,? 

^£=0 V Y^Zo e 2l3mmz ~ pmz2 s(s(i-M),e{z-M)) J 



1=1 ,/^£=0 ^ 

(15) 

Then, the average case performance is determined by the following averaged 
mean square error: 



1 N 

(16) 

Solving these self-consistent equations with respect to m (14) and m (15), 
we evaluate the statistical performance of the present method through the 
quantity a (16) analytically. 

As we have estimated using the Monte Carlo simulation, we estimate how the 
mean square error depends on the hyper-parameter T m for the infinite-range 
version of our model when we set to Q — 8, J s — 1, T s — 1, M — 3.5 (= 
(Q - l)/2),4.5 and J = 1. 

We find from Figs. 6 (a) and (b) that the mean square error takes its minimum 
in the wide range on T m including the Bayes-optimal condition T m = T s (= 1). 
Here, we note that m = m Q (= 3.5) holds under the Bayes-optimal condition, 
T rn = T s for both cases of M — 3.5 and M = 4.5, which is shown in Fig. 
7. From this fact, we might evaluate the gap A between the lowest value of 
the mean square error and the second lowest value obtained at the higher 
temperature than T s as follows. 
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Q 2 E\ =0 



- 1 e 2/3 s mo?-/34 2 

E?= 1 (2e - 2m + ly^t-M 2 
Q 2 E^To 1 e 2 ^ m ^-^ 2 



Q 2 E®=o e 2 P° m °t-P°t 2 



1 

Q 5 



(17) 



For example, for Q = 8, we evaluate the gap as A = (8) 2 = 0.00156 and 
this value agree with the result shown in Fig. 6. From Figs. 6 and 7, we also 
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(b) 

Fig. 6. (a) The mean square error as a function of the parameter T m when Q = 8, 
T s = 1, J s = 1, M = (Q — l)/2 and J = 1, (b) The mean square error as a function 
of the parameter T m when Q = 8, T s = 1, J s = 1, M = 4.5 / (Q - l)/2 and J = 1. 
The value mi for each line caption denotes the initial condition of the magnetization 
m to find the locally stable solution. 

find that the range of T rn in which the mean square error takes the lowest 
value coincides with the range of temperature T m for which the magnetization 
satisfies m(T m ) = m(T s = 1) ± 1 = 3.5 ± 1 as shown in Fig. 7. This robustness 
for the hyper-parameter selecting is one of the desirable properties from the 
view point of the practical use of our approach. 
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Fig. 7. The magnetization m as a function of the parameter T m when Q = 8, T s = 1, 
J s = 1, M = 4.5 ^ (Q - l)/2 and J = 1. 

Moreover, the above evaluations might be helpful for us to deal with the inverse 
halftoning from the halftoned image of the standard image with confidence. 
In fact, we are also confirmed that our method is practically useful from the 
resulting image shown in Fig. 5 (right) having the mean square error a = 
0.002005. 



4 Summary 

In this paper, we investigated the condition to achieve the Bayes-optimal per- 
formance of inverse halftoning by making use of computer simulations and 
analysis of the infinite range model. We were also confirmed that our Bayesian 
approach is useful even for the inverse halftoning from the binary dots obtained 
from standard images, in the wide range on T m including the Bayes-optimal 
condition, T m = T s . We hope that some modifications of the prior distribution 
might make the quality of the inverse halftoning much better. It will be our 
future work. 
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